home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / MATH1 / DIFFUS.PAS < prev    next >
Pascal/Delphi Source File  |  1985-04-03  |  4KB  |  142 lines
  1. { die ln-funktion bei MT+ ist fuer Werte kleiner 1E-5 aeussertst langsam,
  2.   so dass man das Gefuehl bekommt, das dass Programm sich aufhaengt. In-
  3.   folgedessen konnte ich dieses Programm nicht verifizieren. -Juergen }
  4.  
  5. program diffus;    { --> 302 }
  6.  
  7. { Pascal program to perform a linear least-squares fit }
  8. { for the diffusion of Zn and Cu }
  9. { with Gauss-Jordan routine }
  10. { Sperate modules needed:
  11.             GAUSSJ,
  12.             PLOT }
  13.  
  14.  
  15. const    maxr    = 20;        { data prints }
  16.     maxc    = 4;        { polynomial terms }
  17.     r    = 1.987;    { gas constant }
  18.  
  19. type    ary    = array[1..maxr] of real;
  20.     arys    = array[1..maxc] of real;
  21.     ary2    = array[1..maxr,1..maxc] of real;
  22.     ary2s    = array[1..maxc,1..maxc] of real;
  23.  
  24. var
  25.     x,y,y_calc    : ary;
  26.     t,d,resid    : ary;
  27.     coef,sig    : arys;
  28.     nrow,ncol    : integer;
  29.     correl_coef,srs    : real;
  30.  
  31.  
  32. external procedure cls;
  33.  
  34. procedure get_data(var x,y,t,d: ary;
  35.            var nrow: integer);
  36. { get values for nrow and arrays t,d }
  37.  
  38. var    i    : integer;
  39. begin
  40.   nrow:=7;
  41.   t[1]:=600.0;    d[1]:=1.4E-12;
  42.   t[2]:=650.0;    d[2]:=5.5E-12;
  43.   t[3]:=700.0;    d[3]:=1.8E-11;
  44.   t[4]:=750.0;    d[4]:=6.1E-11;
  45.   t[5]:=800.0;    d[5]:=1.6E-10;
  46.   t[6]:=850.0;    d[6]:=4.4E-10;
  47.   t[7]:=900.0;    d[7]:=1.2E-9;
  48.   for i:=1 to nrow do
  49.     begin
  50.       x[i]:=1.0/(t[i]+273.0);
  51.       y[i]:=ln(d[i])
  52.     end
  53. end;        { procedure get data }
  54.  
  55.  
  56. procedure write_data;
  57. { print out the answers }
  58. var    i    : integer;
  59. begin
  60.   writeln;
  61.   writeln;
  62.   writeln('  I      TC        D    DCALC');
  63.   for i:=1 to nrow do
  64.     writeln(i:3,t[i]:8:0,d[i],' ',y_calc[i]);
  65.   writeln; writeln(' Coefficients ');
  66.   writeln(coef[1],' constant term');
  67.   for i:=2 to ncol do
  68.     writeln(coef[i]);        { other terms }
  69.   writeln;
  70.   writeln('D0=',(exp(coef[1])):7:2,' cm sq/sec.');
  71.   writeln('Q =',(-r*coef[2]/1000.0):8:2,'kcal/mole');
  72.   writeln;writeln('SRS= ',srs:7:3)
  73. end;        { write_data }
  74.  
  75. {procedure square(x: ary2;
  76.          y: ary;
  77.          var a: ary2s;
  78.          var g: arys;
  79.      nrow,ncol: integer);}
  80. { matrix multiplication routine }
  81. { a= transpose x times x }
  82. { g= y times x }
  83.  
  84. {$I C:SQUARE.LIB }
  85.  
  86. {external procedure gaussj(var b:    ary2s;
  87.                   y:    arys;
  88.               var coef:    arys;
  89.               ncol:        integer;
  90.               var error:    boolean);
  91. }
  92. {$I GAUSSJ.LIB }
  93.  
  94. procedure linfit(x,        { independant variable }
  95.          y: ary;    { dependent variable }
  96.          var y_calc: ary;    { calculated dep. variable }
  97.          var resid:  ary;    { array of residuals }
  98.          var coef:   arys;    { coefficients }
  99.          var sig:    arys;    { error on coefficients }
  100.          nrow:       integer;    { length of array }
  101.          var ncol:   integer);    { number of terms }
  102.  
  103. { least squares fit to nrow sets of x and y pairs of points }
  104. { Seperate procedures needed:
  105.     SQUARE -> form square coefficient matrix
  106.     GAUSSJ -> Gauss-Jordan elimination }
  107.  
  108. var    xmatr        : ary2;        { data matrix }
  109.     a        : ary2s;    { coefficient matrix }
  110.     g        : arys;        { constant vector }
  111.     error        : boolean;
  112.     i,j,nm        : integer;
  113.     see,a1        : real;
  114.  
  115. begin        { procedure linfit }
  116.   ncol:=2;    { number of terms }
  117.   for i:=1 to nrow do
  118.     begin        { setup matrix }
  119.       xmatr[i,1]:=1.0;    { first column }
  120.       xmatr[i,2]:=x[i]    { second column }
  121.     end;
  122.   square(xmatr,y,a,g,nrow,ncol);
  123.   gaussj(a,g,coef,ncol,error);
  124.   srs:=0.0;
  125.   a1:=exp(coef[1]);
  126.   for i:=1 to nrow do
  127.     begin
  128.       y_calc[i]:=a1*exp(coef[2]*x[i]);
  129.       if y[i]<>0.0 then resid[i]:=y_calc[i]/y[i]-1.0
  130.       else resid[i]:=y[i]/y_calc[i]-1.0;
  131.       srs:=srs+sqr(resid[i])
  132.     end
  133. end;    { linfit }
  134.  
  135.  
  136. begin        { main program }
  137.   cls;
  138.   get_data(x,y,t,d,nrow);
  139.   linfit(x,y,y_calc,resid,coef,sig,nrow,ncol);
  140.   write_data
  141. end.
  142.